Analytical and numerical study of uncorrelated disorder on a honeycomb lattice 
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We consider a tight-binding model on the regular honeycomb lattice with uncorrelated on-site 
disorder. We use two independent methods (recursive Green's function and self-consistent Born 
approximation) to extract the scattering mean free path, the scattering mean free time, the density 
of states and the localization length as a function of the disorder strength. The two methods 
give excellent quantitative agreement for these single-particle properties. Furthermore, a finite-size 
scaling analysis reveals that all localization lengths for different lattice sizes and different energies 
(including the energy at the Dirac points) collapse onto a single curve, in agreement with the one- 
parameter scaling theory of localization. The predictions of the self-consistent theory of localization 
however fail to quantitatively reproduce these numerically-extracted localization lengths. 



I. INTRODUCTION 

Anderson localization (AL) of waves [l| in disor- 
dered media is a ubiquitous phenomenon which has been 
observed both for classical and quantum waves, e.g. 
light 0, 0, acoustics 0, @, water waves @, ultracold 
atoms [3, [H, polaritons [31 and quantum Hall system [Ioj |. 



The scaling theory of localization predicts that a 
three-dimensional (3D) system exhibits a metal-insulator 
transition while ID and 2D systems always display local- 
ization at any finite disorder strength. Approximate an- 
alytical expressions for the localization length in terms 
of the transport mean free path can be derived within 
the framework of the self-consistent theory of localiza- 
tion [l2|, [TH]. Dimension two is in fact the critical di- 
mension for AL and symmetry considerations can play 
an important role. Indeed, while localization is expected 
to take place for spinless time-reversal invariant systems 
(albeit with an exponentially large localization length), 
perturbative renormalization group studies on non-linear 
tr-models suggest that a metal-insulator transition may 
occur in 2D if chirality is present [HI, EH ■ Such a dis- 
ordered system with different chiral classes could be re- 
alized with the honeycomb lattice. The successful iso- 
lation of graphene flakes in 2004 [l6|], and the discovery 
that graphene samples exhibit a finite electronic conduc- 
tivity at half-filling although the density of states (DoS) 
vanishes [13, EH , has thus spurred interest in studying 
electronic transport in graphene in the presence of disor- 
der, see [HI and references therein. 

However, even though graphene is a readily-available 
physical realization of a honeycomb lattice, its properties 
are invariably affected by the combined effects of interac- 
tion, disorder and phonons. The controlled study of dis- 
order alone in graphene sheets is thus difficult, notwith- 
standing the fact that engineering disorder with given 
statistical properties seems out of reach. In that respect, 
ultracold atoms loaded on a graphene-like optical lat- 
tice (20l - [22^ offer an alternative route and have already 
proven their key impact in weak and strong localization 



studies @, l23l - [2"7l | . Furthermore, key transport quanti- 
ties, like the scattering and transport mean free paths or 
the localization length, have already been analyzed for 
speckle optical potentials in the Born approximation in 
Ref. fl3j while engineering disorder with different corre- 
lation properties is possible. 

The aim of this paper is to study transport in a dis- 
ordered honeycomb lattice. We first stick to the simpler 
case of the tight-binding model for the regular honey- 
comb lattice in the presence of uncorrelated on-site disor- 
der characterized by a symmetric box distribution. The 
more interesting (but more complicated) case of corre- 
lated on-site disorder will be the scope of a forthcoming 
publication. The novelty of our work lies in the gen- 
eralization of two known methods, the recursive Green's 
function method and the self-consistent Born approxima- 
tion, (i) to extract single-particle properties, such as the 
scattering mean free time r, the scattering mean free path 
£, the density of state (DoS) v, and the localization length 
£, and (ii) to relate them to experimentally-controlled 
parameters, such as the tunneling amplitude J and the 
disorder strength W. The numerical data obtained from 
these two methods show remarkable agreement and give 
an accurate estimation of these single-particle properties. 
Our results further confirm the one-parameter scaling 
hypothesis [lH for localization but also reveal a quantita- 
tive disparity with thepredictions of the self-consistent 
theory of localization [lj, [l3j . This disparity is yet to be 
understood. 

Currently, there are three pieces of numerical works 
that are technically relevant to ours. Using a transfer ma- 
trix technique, Schrcibcr and Ottomeier |28| have shown 
that the localization lengths for various lattices (square, 
honeycomb and triangular) at the energy band centre 
E — - thus including the Dirac point of the honeycomb 
lattice - and for various disorder strengths obey the same 
scaling laws. Using the same method, Xiong and Xiong 
concluded that all states are localized but found that the 
scaling behavior at the charge neutrality point is different 
from the one at different energies [29] . On the other hand, 
Lherbier et al. considered the time dynamics of a ran- 
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dom phase wave packet using a real space order- TV Kubo 
method [3(j. They subsequently extracted the diffusion 
constant, and hence the scattering mean free path (which 
coincides with the transport mean free path as scattering 
by our <5-correlated potential is isotropic) , from the time 
evolution of the spatial spread of the wave packet. This 
extracted scattering mean free path was then used to 
deduce the semi-classical conductivity through Einstcin- 
Kubo relation [3l( and the localization length through 
the self-consistent theory of localization, a procedure that 
our numerical data question. 

The paper is organized as follows. Sec. [TT] gives the 
essential ingredients of our model and the eigenstruc- 
ture of the disorder-free honeycomb lattice in the tight- 
binding regime. In Sec. IIII1 we introduce the self-energy 
and detail the Born and self-consistent Born approxima- 
tions (SCBA). We further derive analytical expressions 
for the self-energy at some particular energies in the 
weak-disorder limit. In Sec. IIV1 we introduce the recur- 
sive Green's function (RGF) method which is designed 
to compute exact matrix elements of the Green's func- 
tion for large system sizes, with the caveat that actual 
computations can take months on a computer cluster. A 
faster but more restricted variant is the recursive trans- 
fer matrix method. Together with a finite-size analysis, 
these methods allow us to extract the localization length 
of the disordered honeycomb lattice at any given energy. 
In Sec. El we perform a finite-size analysis of the localiza- 
tion lengths. We show that they can be simultaneously 
scaled for all energies, including the charge neutrality 
point. Furthermore, our results indicate that this univer- 
sal curve is valid for all lattice types, all energies within 
the energy band and possibly for all types of uncorrelated 
disorders, which is not surprising from the viewpoint of 
the one-parameter scaling hypothesis. In addition to the 
localization length, we also extract the scattering mean 
free path and the DoS from the recursive Green's function 
method evaluated at complex energies. The extracted 
quantities show remarkable agreement with our results 
from the self-consistent Born approximation. The com- 
parison of the numerically computed localization length 
to the prediction of the self-consistent theory of local- 
ization shows a fair qualitative agreement, but marked 
quantitative differences. We finally conclude in Sec. I VII 
Additional details are given in the Appendices. 



II. MODEL 




FIG. 1. A honeycomb lattice with lattice constant a and 
its diamond-shaped two-point basis cell (dashed line). The 
vectors ci (I — 1,2,3) connect an A-site to its three B-site 
nearest-neighbors. 



scale. Throughout the paper, we assume the diagonal el- 
ements Si to be independent random variables character- 
ized by the same symmetric box probability distribution 



w for N ^ ¥ > 

otherwise. 



(2) 



where W is the disorder strength. The disorder has thus 
zero mean average el = and its two-point correlator is 
then Cij = Si£j 



^frSij, where the bar denotes averag- 



ing over disorder configurations. The disorder being spa- 
tially (^-correlated, scattering is isotropic; the scattering 
and transport mean free paths are consequently equal. 

We now briefly review the eigenstructure of the 
disorder- free Hamiltonian Hq. We refer to Ref. [2(J for 
more details. The regular honeycomb lattice being a tri- 
angular Bravais lattice with a two-site basis cell, it can be 
pictured as two shifted triangular sublattices denoted by 
A and B, see Fig. [T] As a consequence, the coordination 
number of the honeycomb lattice is three. We denote by 
c/ (I = 1, 2, 3) the link vectors connecting a site i E A to 
its three nearest-neighbors ji € B (|c/| = a, a being the 
lattice constant). We next define the structure factor of 
the honeycomb lattice for nearest-neighbor hopping as 



We consider here a tight-binding Hamilton operator 
acting on a regular honeycomb lattice with on-site disor- 
der 

H = H a + V = -J + + $>|i><<|, (1) 



where |i) refers to a Wannicr state localized on site i and 
denotes a sum over nearest neighbors. The hopping 
parameter J is usually positive and it fixes the energy 



/(fe) = |/(fe)| e ^ fc )= J> ifc " 



(3) 



i=i 



For the honeycomb lattice depicted in Fig. [TJ where C\ 
points along the y-axis, we have 



l/(*0| J 



l+4cos z 



-4 cos 



V3k x 
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FIG. 2. The first Brillouin zone B (in grey) of the honey- 
comb lattice shown in Fig. [1] and the associated reciprocal lat- 
tice vectors (arrows). The two empty circles mark the Dirac 
points, around which the energy dispersion is linear. By using 
convenient reciprocal lattice translations, the first Brillouin 
zone can be mapped onto the shown rectangle which is then 
used as the integration space in all integrals in the paper. 



As shown in Ref. 
fined by 



I20f . the eigenstructure of Ho is de- 



(5) 
(6) 



H = ^2 e ks \ks)(ks\, 

k£B.s=± 

£ks =se k = sJ |/(fc)|, 

ifc S ) = -i=[E eifc ' r, i j )- se ^ (fc) E^'Ij')] ( ? ) 

where B is the first Brillouin zone of the honeycomb lat- 
tice (see Fig. [2j, s is the band index (s = +1 for the 
upper, or conduction, band; s = — 1 for the lower, or va- 
lence, band) and N c is the number of Bravais cells of the 
lattice. 

The full spectrum span the energy interval [—3 J, 3 J] 
and band crossing can only occur at zero energy since 
£+(k) = £-(fc) implies that f(k) = 0. This happens at 
the Dirac points K and K' = K where 



K 



( 47T 



\3\/3a 







(8) 



for the honeycomb lattice in Fig. [TJ In the solid-state 
community, the energy at the Dirac points is usually re- 
ferred to as the charge neutrality (energy) point because 
the number of energy states above and below this point 
are equal. As a consequence, when the gate voltage is 
fixed at e±(K) = in a graphene sample, the graphcne 
sheet is charge neutral since the particle and hole states 
are balanced. 

In the rest of the paper, we will be mainly interested 
in the different localization properties of our model for 
energies E lying near the band edges, E sa ±3 J, and 
near the band centre, E w 0. In the first case, k lies 
near the centre of the Brillouin zone (ka <C 1, where 
k = \k\), while in the second case k lies near the Dirac 
points (qa <C 1, where q = k — K, and similarly around 
K>). 

Near the band edges E = ±3 J, the dispersion relation 



is quadratic, thus representing free massive particles 

a 2 k 2 ' 



E ks « 3Js 1 



ka < 1, 



(9) 



the effective mass m s = — s2h 2 /(3Ja 2 ) being negative 
for the upper band and positive for the lower band. 

Near the charge neutrality point E = 0, the disper- 
sion relation is linear, thus representing the celebrated 
"rclativistic" massless Dirac particles propagating with 
a velocity c playing the role of an effective speed of light, 



3a J 



(ga< 1), 



(10) 

(11) 



III. SELF-ENERGY 

An initial Bloch state \ks) propagating in the lattice 
will suffer scattering by the disorder fluctuations and thus 
will be depleted, decaying exponentially over time with 
a time constant which is the scattering mean free time 
T ks . This coherent propagation and decay are described 
by the disorder-averaged Green's function G which obey 
the Dyson's equation 



G(z) =G (z) + G (zMz)G(z), 



(12) 



z being a point in the complex energy plane. Within our 
model, the disorder-free Green's function is given by 



G (z) 



1 



\ks)(ks\ 



Z - H ^ Z - Sks 
ks 



(13) 



The Dyson equation features the self-energy £ which is 
a central quantity given by a perturbative sum of irre- 
ducible diagrams [32(. As the disorder average restores 
the lattice translation symmetries and characteristics, 
one has 

(k's'\G(z)\ks) =G ks (z) <W 5 SS ,, (14) 
(k's'\T,(z)\ks) = Efc s (z) S kk > 5 SS >, (15) 

where the dependence of Y, ks (z) on k is usually smooth. 
Due to particle-hole symmetry, the spectrum of Ho is 
symmetric with respect to E = 0. Since the on- 
site energies are themselves independent symmetrically- 
distributed random variables, it is easy to show that 
(ks\G(z)\ks)* = -(k, — s\G(— z*)\k, — s), where the star 
stands for complex conjugation. In turn, the self-energy 
satisfies 



£ fe , s (z) = T,* k J-z*). 



(16) 



1 These identities no longer hold for a speckle potential since it 
breaks the V — ¥ — V symmetry. 
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Furthermore, because of time-reversal invariance, the 
self-energy is the same at ±fc. This means that it is suf- 
ficient to study the negative energy sector and forward 
propagation (k x > 0). The scattering mean free time, 
defined through 



= - liaS ka (E), 



2r ks (E) 

is independent of the band index 



(17) 



A. Born approximation 

An analytical expression for E& s (z) is generally not 
available and one has to resort to approximations to find 
the self-energy. For weak disorder, the simplest approxi- 
mation is the Born approximation which consists in dis- 
carding all terms of the full diagrammatic perturbative 
expansion in Fig.|3]cxcept the first one. Its lattice matrix 
elements are 

(i\Z BoIa (z)\j) = C iS (i\G (z)\j) = I(z) Sn (18) 
with 

f dk z 

^) = (i\Go( Z )\i) = j B ^ z2 _ j2{mr (19) 

where ft = 8n 2 / (3y/3a 2 ) is the area of the first Brillouin 
zone. For uncorrelatcd on-site disorder, the self-energy 
at the Born approximation is a scalar: 

W 2 

£ fcs (z) PS X BoIn {z) = - ( 20 ) 

The average Green's function then reads 

GWn(z) = Tt ~F TT = G o( z ~ S B orn(^))- 

Z- Mo — ^Born(z) 

(21) 

The expression of I(z) in terms of elliptic integrals [33|, 
HH is given in Appendix lAl 

B. Self-consistent Born approximation (SCBA) 



FIG. 3. The "rainbow" subclass of diagrams retained to 
compute the self-energy in the self-consistent Born approx- 
imation. The double line with arrow denotes an averaged 
Green's function G while a single line with arrow denotes a 
disorder-free lattice Green's function Go- Two vertices (solid 
dots) connected by a dashed line represent a 2-point correlator 
dj — &[e]. The Born approximation consists in computing 
the self-energy with the first "rainbow" diagram only. For un- 
corrected disorder, the connected vertices correspond thus to 
the same site. In this case the self-energy is a scalar operator 
depending only on the energy and the disorder strength. 

enforces < 9 < tt. Eq. (TIB)) then translate into 

1 (-E,W)=j(E,W) (23) 
6(-E, W) = tt - 6(E,W). (24) 

This implies 0(0, W) = tt/2 for any disorder strength W . 
Figs. H] and [5] show the independence of 9 and 7 in the 
SCBA for some particular disorder strengths W, while 
Figs. [5] and [7] show the IT-depcndcnce of 9 and 7 in the 
SCBA at some particular energies E. In the following 
subsection, we investigate SCBA analytically in the weak 
disorder regime. 



C. Weak disorder limit 

In the weak disorder regime W <C J, one expects 7 <C 
1. Several analytical results can then be derived in this 
limit. Some details are exposed in Appendix [Bl 



This approximation scheme builds on the Born approx- 
imation by replacing Go by G in Eq. (fT8|) . It is more 
powerful as it amounts to sum the infinite subclass of 
"rainbow" diagrams given in Fig. [3] It gives the follow- 
ing self-consistent equation: 

W 2 — W 2 
Sscba(^) = -yj (i\G(z)\i) = — I(z - S S cba(^)), 

(22) 

which is easy to solve numerically. 

In the following we will use the paramctrization 
SscbaW — "fJe~ w with 7 positive. For the scatter- 
ing time to be positive, we must have ImEfe s < 0, which 



1. Charge neutrality point 

At the charge neutrality point, we have z = E = 
and the function / in (j2"2"j) is thus evaluated at the di- 
mensionless complex number Z = — Escba(0, W)/ J. 
Furthermore, because of Eq. (fTo) . 9(0, W) = tt/2, and 
^SCBa(0, W) = — 17(0, VF) J is purely imaginary. Equa- 
tions (|2"2")l and (|A3[) have two solutions. One is the trivial 
solution 7(0, W) = while the nontrivial solution solves 

[ dk 1 _ 12J 2 

J B IT |/(fc)| 2 + 7 2 (o, w) ~ ~w 2 ~' (25) 
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E/J 

FIG. 4. (Color online) Using the parametrization 
£(z) = 7Je~ ie (7 > 0, < 9 < n), the figure shows the 
angle 9 as a function of the energy E in the Born approxima- 
tion (black open circles) and in the SCBA for various disorder 
strengths. In the Born approximation, 9 is independent of W . 
At W — 0.1J, SCBA and the Born approximation are essen- 
tially identical. 




E/J 

FIG. 5. (Color online) Using the parametrization £(z) = 

in 

jJe (7 > 0, < 9 < 7r), the figure shows the amplitude 7 
as a function of the energy E in the SCBA for various disorder 
strengths. As W — > 0, the peak at E = J develops into the 
van Hove singularity. 




FIG. 6. (Color online) Using the parametrization £(z) = 
yJe (7 > 0, < 9 < tt), the figure shows the angle 9 as 
a function of the disorder strength W in the SCBA near the 
charge neutrality point, at the van Hove singularity and at 
the band edge. 




FIG. 7. (Color online) Using the parametrization £(2) = 
7Je — (7 > 0, < 9 < 7r), the figure shows the amplitude 
7 as a function of the disorder strength W in the SCBA near 
the charge neutrality point, at the van Hove singularity and 
at the band edge. 



An expansion of the elliptic integrals appearing in the 
first line of Eq. ([A3]) for \Z\ < 1 leads to 



7(0, W0 « 3 exp ( — 



W 2 ' 



(26) 



Shon et al. found essentially the sam e ty pe of result 
7,/ = e c cxpi-Ai/W 2 ), see Eq. (3.21) in The differ- 
ence between their parameter values and ours arises from 
their introduction of the cut-off energy e c . In our treat- 
ment, we use the exact dispersion relation, beyond the 
linear approximation around the Dirac points, making 



the introduction of an artificial cut-off unnecessary. 



2. Van Hove singularities 

At the van Hove singularities, the disorder-free DoS 
diverges. For our honeycomb tight-binding model, this 
occurs at z = E = ±J [2(| and the function / in (|22|) 
is now evaluated at the dimensionless complex number 
Z + = 1 — Y,scba(J, W)/ J. A small-parameter expansion 
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of equations (|22|) and (|A3I) around Z = 1 now leads to 

^2 / ,O.A~r 72. /fi/W r2. 



I6nj 



( b (^) 



1 In In 



^2 



9(J,W) ~ 2( 1+ 3 (ln(^)-lnln(^)) 



(27) 



(28) 



can be calculated more easily [36[. These sections are 
then "glued together" one after one and Dyson's equa- 
tion [33] is repeatedly used to derive the full Green's func- 
tion in terms of the Green's functions of the smaller sec- 
tions. To study localization, we will use a generalized 
version [H, [39[ of the RGF method. This generalized 
version enables us to extract any lattice matrix element 
of the Green's function conveniently and with high nu- 
merical stability. 



One sees that, at lowest order, the self-energy at the 
van Hove singularities is purely imaginary too, 6 w n/2. 



3. Band 



At the disorder-free band edges, z = E = ±3 J. By the 
same token, we evaluate the function I in (|22p at the di- 
mensionless complex number Z + = 3 — Escba(3 J, W) / J. 
A small-parameter expansion around Z — 3 now leads to 



Applying the RGF scheme to our case amounts to con- 
sider a finite quasi-lD lattice strip and to divide it into N 
vertical slices, the two open ends being along the horizon- 
tal direction, see Fig. [5] Denoting by iijv-i the nearest- 
neighbor tight-binding Hamilton operator for a strip with 
(N—l) slices, the Hamilton operator Hm obtained by glu- 
ing an additional slice N can be split into three terms, 



J SCBA 



(3 J, W) 



V3W 2 
' 48ttJ 



ln(ln^ 



, /192V3VJ 2 

m = 

V W 2 

( 192y/3nJ 2 



W 2 



m)). (29) 



Anticipating results displayed and discussed in Para- 
graph IV CI within the SCBA scheme, the DoS van- 
ishes outside a finite energy band, with a square-root 
behavior near the band edge. This SCBA band edge 
is approximately given by the solution of the equation 
E — ReEgcBA^) = 3 J and does not coincide with the 
exact band edge 3 J + W/2 found for the box disorder. 
Note also that the Lifshitz tail between 3 J and 3J+W/2 
is completely missed by the SCBA scheme. 



H n — H 



N-l 



H 



slice , rrhop . rrhop / qn \ 

N "T U N-1,N "T U N,N-V V JU J 

where N is the nearest-neighbor hop operator con- 

necting sites within slice (N — 1) to sites within slice 
N (and vice- versa for H^ > ^_ 1 ). Since no external gauge 
fields are present in our model, we safely consider the hop 
operators to be real from now on. H^ lcc is the nearest- 
neighbor tight-binding Hamilton operator for the isolated 
slice N before it is stacked to the others. It thus includes 
the on-site disorder diagonal term V in Eq. ([TJ. 



IV. RECURSIVE GREEN'S FUNCTION 
METHOD 

The RGF method is based upon the division of the sys- 
tem in smaller sections for which the Green's functions 

I 



Using Dyson's equation, the submatrix G> „ of the 

full retarded Green's function G (Ar) = (E + i0+ - Hn)- 1 
coupling slice I to slice n at energy E can be obtained 
through the following recursion relations, 



G 



(JV) 



G 



{N-l) 
l.n 



M 1 *- 1 ) fjhop 



f~({N) _ „(JV_l) rrhop 



G 



l.N 

(AO 
N.n 



G 



l,N-l 
(N) 



N-1,N N.N> 



H 



hop 



N.N ^N,N 



(N-l) 



/-A JV - 
-1 <~*N-1 



G 



(N) 
jV,n> 



(31a) 
(31b) 
(31c) 



U N,N — \ & 



;fl+ rrslicO rrhop 

lU - H N - H N N . 



r {N-l) 
-1 "JV-l.JV- 



H 



hop 

N-1,N 



(31d) 



r 



with 1 < I < n < (N 
derivation. 



1), see Appendix ICl for a short 



In the following subsections, we explain how we chose 
to slice the honeycomb lattice for the zigzag and arm- 
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rrhop 

a N-l,N "f a N,N-l 



r hop 



Lm 



N-l N 



FIG. 8. Schematic illustration of the RGF method with a 
strip of square lattice of length Ljv and width Lm- The system 
under consideration is constructed by repeated stacking of 
vertical slices of length Lm. The Green's function at each 
step of the stacking process is calculated recursively using 
Dyson's equation. Knowing the Green's function G^" 1 ' = 
(E+iO + —Hn-i)^ 1 for the system with (N— 1) slices, the hop 
operator -ff^ op = #jv°-i N + ff/v°5v-i coupling slices (N — 1) 
and N, and the Hamilton operator Hff' cc of the isolated slice 
N, one can exactly compute the Green's function G^ N ' for 
the whole system of N slices, see Eqs. (|31[) . 



chair geometries. Once the stacking is properly denned, 
each lattice site i = (n, m) will then be labelled by two 
integers. The first one 1 < n < N corresponds to the 
slice it belongs to from left to right, N being the total 
number of such slices. The second one 1 < m < M corre- 
sponds to its position along the slice from bottom to top, 
M being the total number of sites per slice. For both 
geometries we have checked that the matrix elements of 
the Green's function obtained through the recursive algo- 
rithm agree well with those computed by direct inversion 
of the operator (E + i0 + — H N ). 



A. Zigzag configuration 

We first consider the zigzag (ZZ) configuration de- 
picted in Fig. |9l where L vertical zigzag chains of the 
honeycomb lattice are stacked along the horizontal di- 
rection. In the following we will be essentially interested 
in the case of periodic boundary conditions in the vertical 
direction and open boundary conditions in the horizon- 
tal one. This imposes the number of sites in each zigzag 
chain to be an even integer 2M. For each zigzag chain, 
one can define A-typc (resp. B-type) vertical slices con- 
taining only A-sites (resp. B-sites). Each of these vertical 
slices contain M sites and there are N — 2L such slices, 
A- type slices alternating with B-type slices. Lattice sites 
are thus parametrized by (n, m) with 1 < n < N = 2L 
and 1 < m < M, see Fig. [9, The width and length of this 
honeycomb ZZ strip are Lm — s/3Ma and ijv ~ 3iVa/4 
(iV> 1). 

With this slicing choice, the Hamiltonian H^ lcc asso- 
ciated to the isolated slice n is simply a M x M diagonal 




FIG. 9. In the honeycomb zigzag configuration with periodic 
boundary conditions along the vertical direction, L zigzag 
vertical chains, each containing 2M sites, are stacked along 
the horizontal direction (L = 4 and M = 3 in the figure). 
Each zigzag chain defines two vertical slices (dashed vertical 
lines) . These slices are labeled from left to right by the integer 
1 < n < N = 2L. Each slice contains M sites labeled from 
bottom to top by the integer 1 < m < M. The slice gender 
alternates between A-type and B-type. Each site in the lattice 
is uniquely parametrized by the couple of integers (n,m). 



matrix with entries given by the M on-site disorder ele- 
ments {e„. m }. Furthermore, one can easily see that 



rrhop 


= -Jl 


if n 


= [mod. 4] 


rrhop 


= JD 


if n 


= 1 [mod. 4] 


rrhop 

n n,n+l 


= -Jl 


if n 


= 2 [mod.4] 


rrhop 


= JD T 


if n 


= 3 [mod.4] 


rrhop 









(32) 

where the T-superscript means matrix transposition and 
D is the M x M matrix given by 



D 



/l 
1 1 
1 1 



\o o o 



0p\ 






(33) 



1 lJ 



When periodic boundary conditions in the vertical di- 
rection are used (which is our case), p = 1. For open 
boundary conditions in the vertical direction, one would 
have p = 0. 

An astute reader may have noticed that the determi- 
nant of D vanishes for M even and periodic boundary 
conditions. In this case the transfer matrix method [361 ] 
cannot be implemented as it would require the inversion 
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of D or D T . To avoid this pitfall, one can neverthe- 
less always choose M odd. Note however that the RGF 
scheme is perfectly immune to this breakdown and its 
implementation does not suffer any flaw as we have duly 
checked. 



B. Armchair configuration 



n, namely 
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if n is odd, 



if n is even, 



(34) 

where p = 1 for periodic boundary conditions along the 
vertical direction (and M even) and p = for open ones. 




FIG. 10. In the armchair configuration (AC), M horizontal 
zigzag chains, each containing N sites, are stacked vertically 
(M = 4 and N = 8 in the figure). There are thus N vertical 
slices labeled from left to right by the integer 1 < n < N 
(dashed lines). Each vertical slice contains M sites labeled 
from bottom to top by the integer 1 < m < M. If periodic 
boundary conditions are imposed along the vertical direction, 
M must be even. The site gender within a slice alternates 
between A-type and B-type. Each site in the lattice is uniquely 
parametrized by the couple of integers (n, m). 

We now turn to the honeycomb armchair configura- 
tion where M horizontal zigzag chains, each containing N 
sites, are stacked along the vertical direction, see Fig. [TU] 
When periodic boundary conditions are imposed along 
the vertical direction, M must be even. The left ver- 
tical boundary of the lattice is now reminiscent of the 
shape of an armchair. Using the same recipe, we slice 
the lattice with vertical lines. There are now N such 
vertical slices, each containing M sites. The width and 
length of this honeycomb AC strip are Lm = 3Ma/2 and 
w V3Na/2 (N > 1). 

In the armchair configuration, it is easy to see that 
the M x M hop matrices satisfy H n °£ +1 = H n °f ln = 
— Jl, while H^ 1CC = — JX„ + e„ where e„ is a M X 
M diagonal matrix with entries equal to the M on-site 
disorder {e„. m } in slice n. X„ is a M X M sparse matrix 
that couples each site to its nearest neighbors within slice 



V. NUMERICAL RESULTS 
A. Localization length £ 

Localized wave functions are expected to decrease ex- 
ponentially at large distances. We thus compute the lo- 
calization length along the horizontal direction as 

A = - lim -LlnTr(p G$[G<», (35) 

where the bar indicates the average over disorder config- 
urations. Note that it corresponds to the log-averaged 
transmission in quasi- ID systems which has the nice 
property of bein g a dditive and self-averaging when new 
slices are added [40| . 

This finite-size localization length Am depends on the 
energy E and the disorder strength W, but also on the 
lattice configuration (ZZ or AC) and width Lm- In all 
our simulations, wc have used a sufficiently large number 
of randomly-generated disorder configurations for each 
set of parameters (E, W, M and lattice configuration) 
such that the estimated relative error in computing 1 / am 
is less than 0.2%. Furthermore 1/Am is always com- 
puted with lengths Lm greater than 5Aa/- This means 
that larger number of samples are required as Am in- 
creases. To avoid numerical underflow, a rescaling of 
G[ N N is done through Eq. (|3Tbl) every 10 multiplications 



2 Some authors define the localization length Xm through ln|i/>| 
(as we do here) or through In |V'| 2 - There is a factor of 2 between 
these two possible definitions. 
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approximately. RGF equations require one matrix in- 
version per slice. To speed up the computation for the 
AC lattice we switched to the recursive transfer matrix 
method [3r| where a true matrix inversion is carried out 
only after the transfer matrix equation is applied for 10 
slices. 

Based on ideas of the renormalization group [3(| and 
the scaling theory of localization [ll| , it was conjectured 
that all data points in a Xm/Lm versus £/ 'Lm plot should 
collapse onto the same universal curve, 



Am 
Lm 



= F 



Lm 



(36) 



where the infinite-lattice localization length £ depends 
only on energy E and disorder strength W. Our results 
in Fig. QT] fully supports this conjecture for sufficiently 
large L M . 

When Lm 2> £, the system becomes insensitive to the 
vertical boundary conditions and one expects Am ~ £■ 
The scaling function should thus satisfy F(x) « x for 
small x. We chose the honeycomb AC configuration at 
E — 3 J and W = 10 J as the reference data set for 
this limiting behavior, therefore fixing £ = (2.46±0.01)a 
in this case, see black data at the bottom left corner 
in Fig. 1111 At finite system width, we nevertheless ex- 
pect the computed localization length Am to be slightly 
shorter than the true infinite-lattice localization length 
£. This behavior is accounted for by a Taylor expansion 
of the scaling function F(x) = x — ax 2 + 0(x 3 ). From 
our data we infer a = (1.11 ± 0.01). 

Starting from the situation at small x, the full scaling 
function F(x) is then built by finding, for each of our data 
sets, the corresponding £ allowing to patch them on a 
single smooth curve. This procedure has been done with 
our AC data sets obtained in the range E G [0.1J, 3 J] and 
W G [1.2 J, 10 J]. We have also checked that AC data sets 
obtained for E slightly larger than 3 J and various W, as 
well as ZZ data sets obtained at different energies (E = 
0.4J, J, 2.9 J ) and W = 1.6 J, all collapse onto this very 
same curve too. In particular, we numerically confirm 
that at E = 0.4J, E = 2.9J, and W = 1.6J, we get the 
same £ for the ZZ and AC configurations. This might be 
understood as a consequence of the isotropic dispersion 
relation in the two energy ranges, see Eqs. (J9j> and (|llj) . 
On the other hand, at E = J, the extracted £ are slightly 
different for the two honeycomb configurations. 

For the AC configuration, an important finding is that 
data sets obtained at the charge neutrality point (E = 0) 
for W > 2.5 J can all be scaled by F(x). Fig.[T2lshows the 
collapsed data for W = 2.5 J, 4J and 6 J (£ = 500a, 70.9a 
and 16.9a respectively). Performing similar calculations 
for different lattice models (honeycomb, square and tri- 
angular), Schrciber and Ottomeier also found that their 
data obtained at E = and W > 4J could be col- 



lapsed onto a single universal curve [28[. Since our data 
at W — 4 J and 6 J with M = 46 agree with theirs at 
M = 40 (after proper unit conversion), we infer that we 
indeed found the same universal curve. 



While it is difficult to numerically confirm that data at 
E = and W < 2.5 J can be scaled by F(x) (see analysis 
below for the reason of this difficulty), there is no appar- 
ent reason why they should not. Therefore, our results, 
combined with those in Ref. [28| . imply that all curves 
Am / Lm as a function of 1/ Lm can be collapsed onto the 
same universal curve by an appropriate one-parameter 
scaling, independently of the lattice type and disorder 
strength, at least for uncorrelatcd box-distributed on-site 
disorder and energies within (and slightly outside of) the 
energy band of the clean system. This is in marked con- 
trast with the findings of Ref. [29| , which claim a different 
scaling behavior at E = (see below for a possible ex- 
planation of this apparent contradiction) . 

For the clean system near E = 0, the energy shell 
slices the band structure near the Dirac points, defin- 
ing two circles of allowed wavevectors centred on K and 
—K. The initial state \ks) belongs to one of these cir- 
cles, say the one around K. Introducing the deviation 
wavevector q = k K , the circles have radius \q\. In the 
presence of disorder, the two circles are broadened in the 
energy shell and form rings of allowed wavevectors with 
mean radius \q\. Since disorder is uncorrelated in space, 
the two rings are coupled by scattering. However, the 
scattering processes being elastic, starting from the ini- 
tial state q, only these two rings around the Dirac points 
will be populated in the course of time. Loosely speak- 
ing, if scattering does change the direction of q, it cannot 
change its modulus. Around E = 0, the dynamics is thus 
characterized by small wavevectors, i.e. by long distances 
in real space. This is why we expect finite-size effects to 
be larger around E = 0. This can be seen in Fig. [13] 
where the dependence of Am on Lm shows a damped os- 
cillatory behavior when Lm/cl is not large enough. As 
explained below, the very existence of oscillations can be 
traced back to the opening of new scattering channels 
when Lm increases whereas the damping can be traced 
back to disorder broadening. 

Consider the finite-size AC configuration. For pe- 
riodic boundary conditions, the allowed values for the 
wavevector k along Ox and Oy are quantized according 
to k x = n x Ak^ and k y = n y AkM with Afcjy = 27r/Ljv, 
Ak M = 2n/L M , n x G {(),■■ ■ ,N} and n y G {0, ■ • • ,M}. 
In the vicinity of E = 0, the curve at constant E will 
enclose a small circular area around the Dirac point K 
in the Brillouin zone containing few such discrete points. 
When the system size, notably Lm, is increased, more 
points enter this area, which means that more channels 
open for propagation. The period of the oscillations ob- 
served in Fig. [T3] thus corresponds to the change in Lm 
allowing one new channel to open, i.e. ALm = 27r/g, 
where g is the radius of the circle at energy E,. Us- 
ing Eq. (fTTj) , we estimate the period of the oscillations 
to be roughly of the order of ALm = 3iraJ/\E\. How- 
ever, each quantized fc-mode corresponding to energy E 
is broadened by disorder, typically by i~ x . As a conse- 
quence modes cannot be distinguished from each other 
when Afcjv/£ w 1. Near the Dirac point, one has £ = or 
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FIG. 11. (Color online) Single-parameter scaling law for the AC honeycomb lattice with uncorrelated on-site disorder. We 
numerically compute the finite-size localization length Am at various energy 0.1J < E < 3J and disorder strength 1.2J < W < 
10J. The fact that all data sets can be collapsed onto a single log-log curve where £ is the infinite-lattice localization length 
confirms the scaling theory of localization [Tj]. Each disjoint color represents a parameter pair (E, W) as shown in the inset. 



and we find that the oscillations due to the opening of 
new channels are washed out as soon as 

. (37) 

LiM 07TJ 

Using SCBA at W = 1.6 J, this simple estimate gives 
a/L M w 4 x 10" 3 at E = 0.1J and a/L M w 7.4 x 10" 3 
at E = 0.2 J, in perfect agreement with the data collected 
in Fig. [□] 

To obtain a reliable estimate of £, it is necessary to 
go beyond the oscillatory region. This is numerically 
challenging at E = for weak disorder, sec (f2"B")) . For 
example, the SCBA estimate at W = 1.6J now gives 
a/LM ~ 2 x 10 -6 , way out of the capabilities of the RGF 
scheme. We believe this is the reason why a different 
scaling at E = was claimed in [2!|: the sample sizes 
used by the Authors were probably not large enough to 
reach the scaling region. As seen in Fig. Q2J it is easy 
to miss the oscillations at E = for small W. The curve 
there is clearly different from the other ones and cannot 
be collapsed by translation (in log scale) onto the uni- 
versal curve, F(x). This might be the origin of the claim 



in Ref. [2!| that a different scaling is observed at E = 0. 
This is however only a finite-size effect. 

In Fig. [T4] we show the variations of the localization 
length £ as a function of the energy E for different dis- 
order strengths W . In the linear dispersion regime (but 
not too close to E = 0), £ appears to be a constant 
whose magnitude increases as W decreases. When get- 
ting closer to E = 0, £ decreases (see curve at W = 2.5 J) 
and a small dip occurs. A similar feature is predicted 
in Ref. [30j when computing £ from the numerically- 
evaluated semi-classical conductivity. Extrapolating our 
results, we find that £ is finite at E = but increases 
sharply when W decreases. Taking graphene as an ex- 
ample (a = 0.142nm, J w 2.7eV), the localization length 
£ is of the order of 10 6 a w 0.1mm at W = 1.6 J w 4.32eV. 
This sample size is achievable with the current state-of- 
the-art technology. 

A second notable feature is the small dip slightly be- 
low the van Hove singularity at E = J. This is surpris- 
ing as the DoS at the van Hove singularity diverges in 
the absence of disorder, which implies more propagation 
channels. Finally we see that £ decreases as E gets closer 
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FIG. 12. (Color online) Collapsing AC data obtained at E = 
and disorder strengths W > 2.5J onto the scaling function 
F(x) shown in Fig. QT] For W = 2.5J (circles), 4J (squares) 
and 6 J (diamonds), we have £ = 500a, 70.9a and 16.9a re- 
spectively. The largest transverse sizes are given by M = 300, 
70, and 70 for W = 2.5J, 4J and 6J respectively. Data for 
W = 2.5J display a small oscillation around the universal 
curve. 
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FIG. 13. (Color online) Plot of Am /£m as a function of a/LM 
at weak disorder W = 1.6J. The oscillations observed at 
small Lm results from the opening of new scattering chan- 
nels when Lm is gradually increased. In the linear dispersion 
regime, the oscillation period in Lm is crudely approximated 
as inversely proportional to E, with a small correction due to 
disorder. 



to the band edge (E = 3 J) and slightly extends beyond 
it. 

To conclude this subsection, we compare in Fig. ll4l our 
numerical data near E = and E = 3 J at W = 1.4 J with 
the analytical prediction of the self-consistent theory of 
localization [12| (see section |V B 3[) . As can be seen, if 
the general trend is qualitatively satisfactory and even 
semi-quantitatively correct (within a factor 10) near the 
band edge, the prediction is off by more than 3 orders of 
magnitude near the band centre E = 0. 



FIG. 14. (Color online) Infinite-lattice localization length £ 
(in units of a) versus energy E (in units of J) for various 
disorder strengths W. For each set (E,W), £ is extracted 
from our numerically-computed Am, Eq. (J35J) , with the help 
of the scaling function F(x) shown in Fig. 1111 The orange 

solid line gives £ = H\J e 7rre ^ — 1 at W = 1.4J as inferred 
from the self-consistent theory of localization, see Eqs. (|45[) . 



In this analytical prediction, the scattering mean free path I 
is estimated using SCBA and n is defined in Eq. (|45[) . If the 
qualitative behavior is relatively satisfactory, the quantitative 
predictions are off by several orders of magnitude in the linear 
dispersion regime and by about an order of magnitude in the 
quadratic dispersion regime. 



B. Scattering mean free path £ 

1. RGF numerical estimate 

The elastic scattering mean free path I defines the dis- 
tance a particle travels on average without being scat- 
tered. For the negative energy sector (see the reason why 
in the following subsection), we compute the ID averaged 
retarded wave function as 



M 

E 

m,m' — 1 



G 



(N+2n) 
n,n+N 



(38) 



where the retarded submatrix G^ 1 ^^ connecting slices 
n and (n + N) is evaluated at complex energy E + irj 
(E < 0) and disorder strength W, r\ being a small positive 
number. The indices m and m! label the sites within each 
slice respectively. In the following, n is chosen such that 
the distance between slice n and the closest open-end 
boundary is large enough (sec below). 

As shown in Fig.[15l | \I/ jv | 2 falls off exponentially with 
a decay constant t(rj). The scattering mean free path is 
further obtained as i = lim^rj ((i])- In our simulations, 
the number of disorder configurations has been chosen 
such that we have 2 x 10 5 samples near the band edge 
at — 3J (quadratic dispersion regime) and 2 x 10 6 sam- 
ples near the band centre (linear dispersion regime). We 
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FIG. 15. Disorder- averaged wave function tj) N shown in the 
complex plane as the longitudinal size of the lattice, fixed by 
the number of slices N, increases (see arrow for the flow). 
The plot has been made for the honeycomb AC configuration 
with M — 200 sites per slice. The complex energy is E + 
177 = — 2.9J + i0.05J and the disorder strength is W = 1.4J. 
The spiraling feature shows that ~ exp (ik w {rj)LN) X 



exp 



(__ Ln\ 



from which we can extract the decay constant 



£(r/). The scattering mean free path is further obtained as 
lim^o^). 



have imposed a minimum distance of 5£(rj) between slice 
n and its closest open-end boundary, thereby fixing the 
value of n in Eq. ((55)) . This ensures that the exponen- 
tial decay is insensitive to the open-end boundaries. Our 
numerical results are thus equivalent to those that would 
have been obtained with an infinite tube. We would like 
to mention here an important technical remark. When 
summing the submatrix entries in (J38j , serious numerical 
cancellation errors occur when E and J have same signs. 
In our case J > and the numerical evaluation is stable 
only if E < 0. The reason behind this numerical instabil- 
ity is the relative sign between the quasi-Bloch sublatticc 
components of \ks), see Eq. [7] 

The decay constant £{rf) is extracted by performing 
a linear fit on ln^jvl a s a function of L/y We see in 
Fig- IS] that a/£(rj) displays oscillations similar to those 
observed for £ when the transverse size (fixed by the num- 
ber of sites M) increases. Again, the oscillations can be 
accounted for by the opening of new scattering channels. 
By further increasing M, the oscillations damp and fi- 
nally the fluctuations in ajl{r\) reach the level of the er- 
ror bars themselves. This happens when the width Lm 
roughly exceeds 10^(77). Beyond this point, the system 
does not feel anymore the transverse periodic boundary 
condition and the corresponding £(rj) are reliable esti- 
mates of the infinite-lattice case. We show how to extract 
£ from these reliable estimates in Fig. [T7] Note that a 
direct evaluation of £ at r\ = is in principle possible, 
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FIG. 16. Oscillation of l/£{rf) as a function of the trans- 
verse size M. The oscillation period can be accounted for 
by the presence of new scattering channels. The oscillation 
stops when the system does not sense the periodic boundary 
condition, i.e. Lm > W£(r)). 
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FIG. 17. (Color online) Evaluation of the infinite-lattice scat- 
tering mean free path £ at a fixed energy E and disorder 
strength W (E = -0.4J and W = 1.2J in the figure). A 
sequence of reliable estimates a/l(rj) is first obtained by the 
RGF technique for smaller and smaller r\ (black open circles). 
A quadratic fit (solid black line) is then used to interpolate 
the value at r\ = 0, from which £ is deduced. For comparison, 
we give the results obtained with the SCBA method (open 
red squares). For the chosen disorder strength the agreement 
is excellent. 



but is unfortunately affected by huge fluctuations as the 
Green's function is computed on the real axis where its 
poles lie. This direct method proves thus much less effi- 
cient in practice. 
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2. SCBA estimate 

In the same spirit, an SCBA estimate of i can be ob- 
tained by using Gq{E— in (|3"5|) and by considering 
a finite horizontal lattice strip of length L with periodic 
boundary conditions in the transverse direction in the 
limit Lm,L — y oo. In the limit of an infinite AC hon- 
eycomb lattice, a careful but straightforward calculation 
using Eqs. and ([IB"]) leads to 



JV 



dq 

2tt z + J(l + 2cosg) : 



(39) 



where z = E — Y1(E) and q = \/3k x a/2. The self-energy 
S(-E) is computed using Eqs. (|2"2"|) and (|A3|) . The distance 
between the two slices being Ln = V~3aN/2 in the AC 
configuration, we then get the scattering mean free path 
as 



lim 



In I* 



AM 



Ln 



(40) 



The exponential decay of ^n is given by the imaginary 
part of the complex pole Q = Q r + iQi solving 



cos Q = — 



l + z/J 



(41) 



with Qi > 0. The scattering mean free path is then just 
a/£ = 4Qj/v3. Approximate solutions are given in the 
Appendix [B] 

A word of caution is here necessary. Inspection of (|39[) 
in the absence of disorder (S = 0) shows that the allowed 
energy range is restricted to — 3J < E < J. A part of 
the positive energy sector is thus missed and will keep 
missed at weak enough disorder. As a consequence us- 
ing Eq. (|3"9"]l is only well adapted to the negative energy 
sector. The physical reason is that the honeycomb lat- 
tice has a two-point Bravais cell and the prescription (f3"5)) 
amounts to consider a symmetric combination of ampli- 
tudes associated to A-sites and B-sites within a Bravais 
cell of the initial slice. Would one had chosen the anti- 
symmetric combination, then one would have gotten 



A? 



J 



2jt 



dq 
2^ 



JNq 



J(l + 2cosg) ; 



(42) 



which is well adapted to the positive energy sector. As 
a rule of thumb, we thus only use Eqs. (|38|) and (|39| for 
the negative energy sector E = — \E\ and then resort to 
(fl~6]) whenever necessary. 

Once the scattering mean free path £(E) and the scat- 
tering mean free time t(E) have been calculated, one can 
compute the ratio v(E) = £(E)/t(E). From gl]), we get 



v(E) 2 sinhQ, 



(43) 



where c is the Dirac fcrmions speed and where < Q r < 
27r/3 for propagation from left to right. At sufficiently 



weak disorder, one expects Qi <C 1, and thus v(E)/c 
(2 sin<3r)/\/3, where Q r solves 



~(Qr) = - J(l + 2 cos Q r ) w E - ReE(J5). 



(44) 



In the weak disorder regime, the real part of self-energy 
is generally small compared to the value of E and can be 
usually discarded. 

Returning to fully dimensioned quantities, we thus find 
the usual result that v(E) = \dk T s(k x )\/h is the group 
velocity (here along Ox) when disorder is sufficiently 
weak. Quantitative comparison between the SCBA and 
RGF results are shown in Figs. [TH] and [TU The two 
methods show remarkable agreement. We note however 
that the SCBA underestimates £ in the quadratic disper- 
sion regime but overestimates I in the linear dispersion 
regime. The deviation of the Born approximation from 
the true (RGF) £ might be understood from two perspec- 
tives: first as a consequence of the smoothing of the DoS 
v{E) by disorder (see Fig. [^JJ), second as a consequence 
of the enhancement of the group velocity (see Fig. |2"T1 
and |2"2"| . Indeed, from the relation £ = vt (v is the group 
velocity), we know that l/£ is directly proportional to 
ImS(_E), which is in turn directly proportional to the 
DoS. As disorder is increased, the sharp variations of the 
DoS at the band edges (E = ±3 J), at the van Hove sin- 
gularities (E = ±J), and at the charge neutrality point 
(E = 0) are smoothed. Since the DoS area must be con- 
served (particle number is conserved), this smoothing is 
accompanied by a redistribution of states over energies 
and by an increase of the DoS at some energies. For ex- 
ample, since the van Hove singularity peaks in the DoS 
are decreased, there is a corresponding increase of the 
DoS in the linear dispersion regime. Similarly, as the 
band edges are smoothed, the DoS at \E\ < 3 J decreases 
but increases for \E\ > 3 J. Therefore, the departure 
of the actual l/£ from the Born value as disorder is in- 
creased, and whether the Born value underestimates or 
overestimates the actual l/£ in Fig. ITS1 and Fig. [TH1 is a 
simple consequence of the smoothing effect of the DoS in 
different energy regimes. 



3. Comparison with the self-consistent theory of localization 

The self-consistent theory of localization (SCTL) [l2| 
provides a self-consistent recipe which extends the weak- 
disorder diagrammatic approach to the regime of Ander- 
son localization. The diagrammatic approach describes 
perturbatively the weak localization corrections to clas- 
sical transport due to interference along closed loops. As 
this correction diverges for infinite system size, the self- 
consistent theory of localization aims at computing the 
minimum size of the system beyond which the correc- 
tion is strong enough to stop the diffusive transport and 
identifies this length scale with the localization length. 

For isotropic scattering and an isotropic dispersion re- 
lation, SCTL establishes a simple link between the local- 
ization length £, the scattering mean free path £ and the 
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FIG. 18. (Color online) Inverse of the scattering mean free 
path £ (in units of the inverse of lattice constant a) extracted 
from the recursive Green's function (RGF) method and from 
the self-consistent Born approximation (SCBA) as a func- 
tion of the square of the disorder strength W in units of 
the hopping energy J. To compare results obtained for the 
honeycomb and square lattices, (W/ J) 2 is renormalized by a 
factor proportional to the corresponding densities of states, 
Vho = 0.14 (honeycomb) and i/ sq = 0.081 (square), in the ab- 
sence of disorder as defined by Eq. (|47|l. For both lattices, the 
(negative) energy E is chosen near the band edge where the 
dispersion is quadratic. The inset identifies the different lat- 
tices, configurations, energies and calculation methods. The 
overlap between the results of the armchair (AC) honeycomb 
and square (SQ) lattices show that the average propagation 
near the band edge is independent of the lattice type and 
solely determined by the quadratic nature of the dispersion 
relation. Comparison with the Born approximation is shown. 
We also plot the scattering mean free path corresponding to 
the numerically observed localization length (calculated with 
the Recursive Green Function method), assuming that the 
two quantities are connected by eq. (|45p derived from the 
Self-Consistent Theory of Localization. This shows that the 
predictions of the Self-Consistent Theory of Localization are 
qualitatively correct, especially for weak disorder, but that 
large quantitative deviations are observed at large disorder. 



wavcvcctor k, 



£ = 2£ \/e™ £ - 1 



(45) 



where k = |fe| when the energy is chosen near the band 
edges and k = \q\ = \k — K\ when the energy is chosen 
near the charge neutrality point. SCLT being valid when 
tit 3> 1 , we further get the approximation 

^ « ln(7r</4) - m(ln(7T</4)). (46) 

We show this SCTL prediction in Fig. HU and Fig. [H 
£ being evaluated by the RGF method. Much to our 
surprise, the agreement between our numerical data and 
the SCTL prediction proves very poor. It should be re- 
minded however that SCLT is supposed to be valid when 
k£ 3> 1, that is in the weak disorder limit where a/£ is 
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FIG. 19. (Color online) Inverse of the scattering mean free 
path I (in unit of the inverse of the lattice constant o) as a 
function of the square of the disorder strength W in units 
of the hopping energy J for the AC honeycomb lattice. The 
(negative) energy E = — 0.4J is chosen in the linear disper- 
sion regime. Connected filled black circles: RGF results. Blue 
dashed line: SCBA results. Black solid line: Born approxi- 
mation (BA). Connected filled red squares: scattering mean 
free path estimated from the localization length, assuming 
that they are connected by eq. (j45]) derived from the Self- 
Consistent Theory of Localization. While the SCBA predic- 
tion agrees well with the RGF numerical computation of the 
mean free path, the predictions of the Self-Consistent Theory 
of Localization are quantitatively off. 

simply proportional to W 2 . Fig. [TBI (and to a lesser extent 
Fig. [19] too) indeed shows that the SCLT estimate tends 
to be a linear function of W 2 . In fact all predictions, 
SCBA, RGF, SCTL and Born approximation agree well 
when W -C J ■ It would be desirable to have numerical 
results at smaller disorder strengths. Unfortunately, the 
localization length is too large to be measurable. In any 
case, our numerical results show that higher orders in W 
are completely different for the true (RGF) a/ 1 and the 
SCLT prediction. Furthermore, the trend in the linear 
dispersion regime is completely different. These results 
are yet to be understood. 



C. Density of states 



by 



The disorder-averaged DoS per lattice site is defined 



{ei j, w/j) = E S ( E I J - V J )> 



(47) 



where iV c is the total number of Bravais cells and where 
is the eigenvalue of the Hamilton operator H in 
Eq. ([T]). This definition is related to the diagonal ele- 
ments of the Green's function by 

v{E/J, W/J) =-- lim lm(j\G(E + i v )\j), (48) 

7T »j->0+ 
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FIG. 20. (Color online) Density of states (DoS) v(E/J, W/J), 
Eq. (|48[) as a function of (negative) energy E (in units of the 
hopping energy J) for several disorder strengths W. The 
number of sites in the transverse direction is M = 200. The 
DoS is increasingly smoothed by disorder as W is increased. 
Because of particle number conservation, the area under the 
curve is also a conserved quantity and states are simply re- 
distributed over a broader energy range. For example, as the 
smoothed van Hove singularity peak at E = — J decreases, 
the corresponding states are redistributed in the wings as in- 
dicated by the arrows. This translates into an increased DoS 
near the charge neutrality point E = 0. Similarly, states at 
the band edge at E — — 3J are partly redistributed outside 
the energy band of the clean system (E < —3 J), resulting 
in a decrease of the DoS in the quadratic dispersion regime 
E > —3 J. The upper inset shows the variations of the DoS in 
the quadratic dispersion regime while the lower inset shows 
the variations of the DoS in the linear dispersion regime near 
the charge neutrality point. 



where j labels an arbitrary lattice site in the infinitely- 
large lattice. Like what we did to extract £ with the RGF 
method, we compute G(E + irj) for r]/J e [0.01,0.07] and 
then perform a quadratic fit in rj to extract the limit 
i] — >• + . For each value of 77, the longitudinal length 
Ljv and the transverse width Lm are both chosen to be 
greater than W£(rf). The value of (j\G(E + irj)\j) is then 
estimated by considering any site j at a minimum dis- 
tance of 5£(rj) from the two ends of the tube. 

The comparison between the DoS calculated with the 
RGF and SCBA methods is shown in Fig. [23] The two 
methods agree remarkably well except near the band 
edges (E = ±3/) and at the charge neutrality point 
(E = 0). The breakdown of SCBA is probably due to 
the fact that q£ <C 1. For moderate disorder strength 
(W = 1.4J), the RGF results display Lifshitz tails near 
the band edges while the SCBA results show a square- 
root cut-off. 

The DoS at E = as a function of W is shown in 
Fig. [Ml The large deviation between the RGF and SCBA 
results is expected since the latter is strictly not valid in 
the range of disorder strengths shown. Note that Shon 
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FIG. 21. Variation of the group velocity v = Ijr with respect 
to the square of the disorder strength W in units of of the 
hopping energy J within the SCBA. The (negative) energy 
E = — 0.4J is chosen near the linear dispersion regime. The 
group velocity of the clean system (chosen along Ox) is vo = 
1.1c, where c is the massless Dirac particles velocity. As one 
can see, the group velocity near the charge neutrality point 
is only slightly enhanced by disorder. The reason why SCBA 
overestimates £ near the charge neutrality point compared to 
the RGF calculation is thus essentially a consequence of the 
smoothing of the density of states by disorder, see Fig. 1201 
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FIG. 22. Variation of the group velocity v — £/t with respect 
to the square of the disorder strength W in units of of the 
hopping energy J within the SCBA. The (negative) energy 
E — — 2.9J is chosen in the quadratic dispersion regime. The 
group velocity of the clean system (chosen along Ox) is given 
by do = 0.36c, where c is the massless Dirac particles veloc- 
ity. As one can see, the group velocity near the band edge is 
significantly enhanced by disorder. This is the main reason 
why SCBA underestimates £ near the band edge compared to 
the RGF calculation. 
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FIG. 23. (Color online) Density of states v(E/J, W/J) com- 
puted by the recursive Green's function method and by 
the self-consistent Born approximation at two weak disorder 
strengths W. The number of sites in the transverse direction 
is M = 200. Since v(E/J, W/J) is an even function of E, wc 
have plotted the RGF (symbols) and SCBA (continuous line) 
predictions at W = 1.4 J in the negative energy sector. The 
corresponding predictions at W — 0.1J have been plotted in 
the positive energy sector. The agreement is excellent. 




W/J 

FIG. 24. (Color online) Density of states v(E/J, W/J) com- 
puted by the recursive Green's function method (full black 
circles) and by the self-consistent Born approximation (rod 
squares) as a function of the disorder strength W (in units of 
the hopping energy J) at the charge neutrality point E = 0. 
SCBA underestimates v(E / J,W / J) by more than a factor 2 
when disorder is large enough. This is the main reason why 
SCBA overestimates £ (see text). 



and Ando [35| obtained slightly different SCBA results 
for the DoS at the Dirac points probably because they 
used a strictly linear dispersion relation, hence discard- 
ing trigonal warping [ill ] . As already witnessed by the re- 
sults on I, the SCBA underestimates v(E), and hence the 
self-energy 12(E), in the linear dispersion regime. How- 



ever, although quantitatively incorrect, the SCBA qual- 
itatively captures an essential property of the system, 
namely the very fast decrease of the DoS when W goes 
to zero, essentially like exp(— A/W 2 ). 



VI. CONCLUSION 

In this paper, we have studied coherent transport 
in the honeycomb lattice subjected to the effect of a 
spatially-uncorrelated on-site disorder with symmetric 
box-like distribution. We have used the recursive Green's 
function method to reliably extract the scattering mean 
free path t and the density of states. We have compared 
these quantities to the analytic predictions of the self- 
consistent Born approximation and found good agree- 
ment at weak disorder. We have also used the recur- 
sive Green's function method to extract the localization 
lengths for different transverse sizes. We have shown 
that all of these finite-size localization lengths can be 
collapsed onto a single curve. We have checked that this 
curve is universal as it applies equally well to the square 
and honeycomb lattices at any energy. In particular, it 
applies to the honeycomb lattice at the charge neutral- 
ity point, at the van Hove singularities, and at the band 
edges. These findings validate the one-parameter scaling 
hypothesis which is thus not restricted to particles with 
either quadratic dispersion or linear dispersion relations. 
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Appendix A: Diagonal elements of the clean lattice 
Green's function 



The diagonal elements of the disorder-free lattice 
Green's function Go turn out to be independent of the 
site i 



I(z) = (i\G (z)\i) 



dk 



n z 2 -j 2 |/(fe)p 



(Al) 



with n = 8tt 2 / (3y/3a 2 ) the area of the Brillouin zone. 
We use the rescaling s/3k x a/2 — s- a and 3k y a/2 —> /3. 
Defining Z = z/J and I(z) = H(Z)/J, little algebra 
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FIG. 25. The boundary Imp 2 = (open circles) divides the 
first quadrant of the complex- Z plane in two regions. The left 
region satisfies KeZ < 1 and Im/i 2 < 0. Because of the analyt- 
ical continuation of the elliptic integral across the boundary 
Im/i 2 = 0, different analytic expressions have to be used to 
compute the rescaled diagonal element of the Green's function 
H(Z) in these two regions, see Eq. (|A3[) . 



where 

g(z) 

M 2 = 



7r(Z-l)i(Z + 3)5 



16Z 



5o(r) = 



pI 



7r(r 2 + i)f (r 2 + 9)3 
16- (y(r 2 + 9)(r 2 + i)- (r 2 + i))' 
4(r 2 + i)5(r 2 + 9)i 



and 



K(p 2 ) = 



dO 



o \J\ - p 2 sin 2 9 



(A4a) 
(A4b) 
(A4c) 

(A4d) 
(A5) 



is the complete elliptic integral of the first kind [42| ■ Fig- 
ure [25] gives the boundary Imp 2 = in the first quadrant 
of the complex- Z plane. 



Appendix B: Self-consistent Born approximation 



1. Van Hove singularities 



gives 



H(Z) = 



* dad/3 



Z 2 — (1 + 4 cos 2 a + 4 cos a cos / 



(A2) 

where Z can be any point in the complex plane except 
the real interval [—3, 3]. 



The idea is now to compute (|A2[) for real Z outside this 
interval and do an analytic continuation in the complex 
plane. In fact, direct inspection shows that KeH is even 
in HeZ and in ImZ, while ImH is odd in HeZ and in 
ImZ. This means that it is sufficient to compute H(Z) 
is the first quadrant (ReZ > 0, ImZ > 0) of the complex 
plane and use these parity properties to infer H(Z) in the 
other quadrants. By demanding I(Z) to vary smoothly 
as Z varies in the complex plane and noticing how p 2 
(see below) crosses the branch cut of the elliptic integral, 

we get m, mi 



H(Z) = 



( iTg {T)K(p 2 ) forZ = ir,r>0 



Zg(Z)(K(p 2 )+2iK(l-p 2 ) 



for Imp 2 < and < RcZ < 1 
Zg{Z)K(p 2 ) otherwise, 



(A3) 



With the parametrization S = "/Jc~ lS , one expects 
7 <C 1 in the weak disorder regime. A small-parameter 
expansion in Eqs. (f2"2")l and (|A3|) . gives at lowest order 
6 w tt/2 and 



A -lJ A - 

7 \7 



64ttJ 2 
~W 2 ~ 



Thus, 



7 



W 2 f6AnJ 2 



16ttJ 2 V ^ 2 



(Bl) 



(B2) 



in terms of the Lambert W -function 0(a) defined for any 
complex number a through the identity 

a = n(a)c Q ( a l (B3) 
At weak disorder, the asymptotic form for |a| 3> 1 [43[ 

' In In a s 



Sl(a) =lna-lnlna + 



In, 



(B4) 



yields Eq. (|27|). 

The self-energy in fact obeys the more general relation 



4J 



^SCBA 



i7(e- i2 -/ 3 (p + m)), 



(B5) 



where p = (1 — i47r/3) and u = 64irJ 2 /W 2 . Assuming 
|p | <C u, we get from the latter equation 6 = tt/2 + 56 
with 



88: 



6(ln(^)-lnln(^)) 



(B6) 
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2. Band edges 

We consider the dimensionless "detuning" S = 3 — E/J 
from the band edge at 3 J. We consider Z = 6 + E/ J 1 
as a small parameter. Similar to what has been done in 
the previous Section, we find that Z approximately obeys 
the equation 



compared to the energy E itself, at least for the disor- 
der strengths considered. This is the case at Ej J = —0.4. 



Appendix C: Derivation of the generalized version of 
the recursive Green's function method 



Z = P(ln(12/Z)-m 



(B7) 



where /3 = Using the Lambert H^-function, we 

find 



-12e 5 /^ 
E m /3Jfi( ) — 5 J. 

P 



(B8) 



The large- argument expansion of the Lambert W- 
function yields 



E « pJ ( ln (^f ) - i7r - m 



]n(j) -m + S/p 



(B9) 



3. Scattering mean free path 

For E < 0, we solve Eq. (|41[) to obtain the following 
approximate solutions: 



4 



1(E) 



E\ < J, 
5 = — J, 



(BIO) 



3 + S/J< 1. 



We note that the factor W 1 — |Re(Z) is needed when 
the energy is sufficiently far from the charge neutrality 
point (where deviations from the linear dispersion regime 
come into play) but not too near the van Hove singu- 
larity (where a different approximation applies). Indeed 
the real part of the self-energy is then no longer small 



Starting with Eq. (|3Tfl) , the Born series for the Green's 
function Gn reads 



Gn = Gn-i + Gn(H 



hop 
N 



rjslicc \ s~1 

n N )Un- 



(Cla) 
(Clb) 



where = H^ 1N + H*°* N _ V In Eq. pi]) , the 

term Gn-iH^^Gn vanishes since H N 1CC couples only 
sites within slice N and one immediately gets Eq. (|31a[) 
after sandwiching with bras and kets of sites not exceed- 
ing the (N— l)th slice. Furthermore, since, by definition, 
Gjv-i does not couple to sites in slice N, we immedi- 
ately recover Eq. (|31b[) by setting n — N in Eq. (|31ap . 
Eq. picp can be obtained from Eq. (|Clb[) through a sim- 
ilar procedure. 

To obtain Eq. (|31d[) . we need to go back to the defi- 
nition of the Green's function (E — Hn)Gn = 1, which 
implies 



(E — ffjv-i — H 



hop 
N 



H 



slice 
N 



)G N = 1. 



(C2) 



Sandwiching Eq. (|C2[) with bras and kets of all sites 
within slice N, we get 



Pr (W) _ rrhop n (N) 
- C/0 JV,JV ■"iV,iV-l u 'iV-l,JV 



rrslice sy(,N) -n 



(C3) 



As Eq. (f3Tb|) yields 



N—1,N " Ur iV-l,JV-l- rI Ar-l,Ar (J JV,JVJ V^^i 

substitution into Eq. (|C3[) immediately gives Eq. (I31dl) . 



-.(JV-l) 



rhop 
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